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-^ Abstract We derive a compatible discretization method that relies heavily on the 

underlying geometric structure, and obeys the topological sequences and commut- 
ing properties that are constructed. As a sample problem we consider the vorticity- 
velocity-pressure formulation of the Stokes problem. We motivate the choice for 
a mixed variational formulation based on both geometric as well as physical ar- 
guments. Numerical tests confirm the theoretical results that we obtain a pointwise 
divergence-free solution for the Stokes problem and that the method obtains optimal 

i-Q convergence rates. 
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1 Introduction 

> 

qq As sample problem we consider the Stokes flow problem in its vorticity-velocity- 

y—i pressure formulation, 

o 

r ■ co-curlu = inf2, (la) 

curla) + gradp = f in Q, (lb) 

^ divu = in 12. (lc) 

• • In this article we consider prescribed velocity boundary conditions, u = on d£2, 

. J^ but the method holds for all admissible types of boundary conditions, see (8). 
S^ Despite the simple appearance of Stokes flow model, there exists a large number 

5— i of numerical methods to simulate Stokes flow. They all reduce to two classes, that is, 
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either circumventing the LBB stability condition, like stabilized methods, e.g. Q, 
or satisfying this condition, as in compatible or mixed methods, e.g. |4|. The last re- 
quires the construction of dedicated discrete vector spaces. Best known are the curl 
conforming Nedelec and divergence conforming Raviart-Thomas spaces. Here, we 
consider a subclass of compatible methods, i.e. mimetic methods. Mimetic methods 
do not solely search for appropriate vector spaces, but aim to mimic structures and 
symmetries of the continuous problem, see [2, 3i [Tol[TTIl . As a consequence of this 
mimicking, mimetic methods automatically preserve most of the physical and math- 
ematical structures of the continuous formulation, among others the LBB condition 
and, most important, a pointwise divergence-free solution, (HE!. 

At the heart of the mimetic method there are the well-known integral theorems of 
Newton-Leibniz, Stokes and Gauss, which couple the operators grad, curl and div, 
to the action of the boundary operator on a manifold. Therefore, obeying geometry 
and orientation will result in satisfying exactly the mentioned theorems, and conse- 
quently performing the vector operators exactly in a finite dimensional setting. In 3D 
we distinguish between four types of sub-manifolds, that is, points, lines, surfaces 
and volumes, and two types of orientation, namely, outer- and inner-orientation. 
Examples of sub-manifolds are shown in Figure [T] together with the action of the 
boundary operator. 
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Fig. 1 The four geometric objects possible in R 3 , point, line, surface and volume, with outer- 
(above) and inner- (below) orientation. The boundary operator, d, maps ^-dimensional objects to 
(k— 1) -dimensional objects. 



By creating a quadrilateral or hexahedral mesh, we divide the physical domain 
in a large number of these geometric objects, and to each geometric object we as- 
sociate a discrete unknown. This implies that these discrete unknowns are integral 
quantities. Since the three earlier mentioned theorems are integral equations, it fol- 
lows for example that taking a divergence in a volume is equivalent to taking the 
sum of the integral quantities associated to the surrounding surface elements, i.e. 
the fluxes. So using integral quantities as degrees of freedom to perform a grad, 
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curl or div, is equivalent to taking the sum of the degrees of freedom located at its 
boundary. 

These relations are of purely topological nature. They form a topological se- 
quence or complex. This sequence is fundamental. It has a direct connection with 
the complexes that are related to the physical domain, the computational domain, 
the physical problem and the discretization. 

Although the original work, [9, 10], was presented in terms of differential geom- 
etry and algebraic topology, here we will use vector calculus because it is the more 
common mathematical language. Nevertheless, we will put emphasis on the dis- 
tinction between topology and metric, on complexes and on commuting diagrams, 
which drives the former two languages. 

We make use of spectral element interpolation functions as basis functions. In 
the past nodal spectral elements were mostly used in combination with Galerkin 
projection (GSEM). The GSEM satisfies the LBB condition by lowering the poly- 
nomial degree of the pressure by two with respect to the velocity. This results in a 
method that is only weakly divergence-free, meaning that the divergence of the ve- 
locity field only convergence to zero with mesh refinement. The present study uses 
mimetic spectral element interpolation or basis functions, [101- The mixed mimetic 
spectral element method (MMSEM) satisfies the LBB condition and gives a point- 
wise divergence-free solution for all mesh sizes. 



2 Can we really discretize exactly? 

Since the Stokes flow model ([TJ) should hold on a certain physical domain, we will 
include geometry by means of integration. In that case we can relate every physical 
quantity to a geometric object. Starting with the incompressibility constraint ( fTc) 
we have due to Gauss' divergence theorem, 

divudV = / undS = 0, 

V JdV 

and using Stokes' circulation theorem the relation (Ta\ can be written as 



0)xndS= /curluxndS= / utd/. 

S JS JdS 

From the first relation it follows that divu is associated to volumes. The association 
to a geometric object for velocity u is less clear. In fact it can be associated to 
two different types of geometric objects. A representation of velocity compatible 
with the incompressibility constraint is given in terms of the velocity flux, u • n, 
through a surface that bounds the volume, while in the circulation relation velocity, 
u ■ t, is represented along a line that bounds the surface. We will call the velocity 
vector through a surface outer-oriented and the velocity along a line segment inner- 
oriented. A similar distinction can be made for vorticity, see |9|. 
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The last equation to be considered is ( fib} . This equation shows that classical 
Newton-Leibniz, Stokes circulation and Gauss divergence theorems tell only half 
the story. From the perspective of the classical Newton-Leibniz theorem, the gradi- 
ent acting on the pressure relates line values to their corresponding end point, while 
the Stokes circulation theorem shows that the curl acting on the vorticity vector re- 
lates surface values to the line segment enclosing it. So how does this fit into one 
equation? In fact, from a geometric perspective, there exists two gradients, two curls 
and two divergence operators. One of each is related to the mentioned integral the- 
orems as explained above. The others are their formal adjoint operators. Let grad, 
curl and div be the original differential operators associated to the mentioned inte- 
gral theorems, then the formal Hilbert adjoint operators grad*, curl* and div* are 
defined as, 

(a, -grad* b) n := (diva,&) fl , (a,curi*b) fl := (curia, b)^, (a,—drv*b) Q :— (grada,b) 

From a geometric interpretation, the adjoint operators detours via the opposite type 
of orientation. Where div relates a vector quantity associated to surfaces to a scalar 
quantity associated to a volume enclosed by these surfaces. Its adjoint operator, 
grad*, relates a scalar quantity associated with a volume to a vector quantity associ- 
ated with its surrounding surfaces. This is illustrated in Figure[2] Following Figure[2j 
the adjoint operator grad* consists of three consecutive steps: First, switch from an 
outer oriented scalar associated to volumes to an inner oriented scalar associated to 
points, then take the derivative and finally switch from an inner oriented vector as- 
sociated to lines to an outer oriented vector associated to surfaces. In a similar way 
we can describe the derivatives curl* and div* . 



Outer Orientation 
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Fig. 2 Geometric interpretation of the action of the boundary operators, vector differential opera- 
tors and their formal Hilbert adjoint operators. 
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Since the horizontal relations are purely topological and the vertical relations 
purely metric, the operators grad, curl and div are purely topological operators, 
while grad*, curl* and div* are metric. This makes them much harder to discretize. 

Now (JTbJ could then either be associated to an inner-oriented line segment by 
rewriting it as 

curl* CO + grad/} = f, 

or be associated to an outer-oriented surface by rewriting it as 

curl co + grad*/? = f. 

Without geometric considerations we could never make a distinction between grad, 
curl and div and their associated Hilbert adjoints div*, curl* and grad*. 

Since our focus is on obtaining a pointwise divergence-free discretization, we 
decide to use the expression where the equations are associated to outer-oriented 
geometric objects, 

CO - curl* u = in Q , (2a) 

curio + grad*/? = f inf2, (2b) 

div u = in Q , (2c) 

where the first equation is associated to outer-oriented line segments, the second to 
outer-oriented surfaces and the third to outer-oriented volumes. 



3 Complexes 

Figure |2]re veals already a number of sequence or complex structures. Starting from 
geometry, we consider points, P, lines, L, surfaces, S, and volumes, V. They possess 
a sequence in combination with the boundary operator, d . The boundary of a volume 
is a surface, the boundary of a surface is a line and the boundary of a line are its two 
end points. This results in the following complex, 

O^pAlAsAv. (3) 

An important property of the complex is that if we apply the boundary operator 
twice, we always find an empty set, e.g. if S = dV, then dS = 0. As follows directly 
from the previously mentioned integral theorems, it follows, as a consequence of 
dd = 0, that curl grad = and div curl = 0. The derivatives themselves also form a 
complex. In a Hilbert setting this becomes, 

H\Q) ^H{Q,cml) ^>H(Q,diw) ^L 2 (Q), (4) 

and using the Hilbert adjoint relations we also obtain the adjoint complex with prop- 
erties curl*grad* = and div*curl* = 0, 
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L\ (fl ) ^ H (Q , div* )^H (n, curl* ) ^ #,} (12 ) (5) 

In the Hilbert setting, the variables of the Stokes problem are in the follow- 
ing spaces, a> e H (Q, curl) nH {£2, div*), u e H(Q, div) n //n (X2 , curl* ) and /? e 
L 2 (f2 ) n Hq (£2 , grad* ) . It is hard, if even possible at all, to find discrete vector spaces 
that are subsets of these function spaces and simultaneously satisfy the complex 
properties. Instead, the Stokes problem can be cast into an equivalent variational 
or mixed formulation where we make use of the Hilbert adjoint properties. This 
simplifies the function spaces of the flow variables. The mixed formulation reads; 



Find (o),u 
for all (d, 


v,tf) 


: {H(Q , curl) x H(Q , div) x L 2 ( 
€ {H(Q , curl) x H{Q , div) x L 


;i2)}wifhfe#(I2,div) 
2 (f2)}, such that, 


given, 






{a,(o) Q - 


-(curlff^)^ 


= 0, 


(6a) 






(v,curl(o) n 


.-(divv,p) n 


= M) fl . 


(6b) 








(^,divu) fl 


= 0. 


(6c) 



With the formulation and corresponding function spaces, we are able to construct 
compatible discrete vector spaces. Note that we now completely avoid the metric 
dependent derivatives grad* and curl* , and their corresponding complex. 

4 Discretization of Stokes problem 

Degrees of freedom. In many numerical methods, especially in finite difference and 
finite element methods, the discrete coefficients are point values. In the proposed 
mimetic structure, the discrete unknowns represent integral values on ^-dimensional 
submanifolds, ranging from points to volumes, so < k < 3. These ^-dimensional 
submanifolds are oriented, constitute the computational domain and span the phys- 
ical domain. The concept of orientation shown in Figures [T] and [2] gave rise to the 
boundary operator, <3, which can be represented by connectivities consisting only of 
-1, and 1, see also @. 

The space of degrees of freedom are given by &, Jz? , 5^ and 1^ , These spaces 
form a duality pairing with the geometric spaces P, L, S and V. The degrees of 
freedom are integral values, i.e. 



w-tdl G_Sf, u-ndSey, pdV e-T. (7) 

; Js Jv 

By the definition of the degrees of freedom spaces and the previously mentioned 
integral theorems, we can define the formal adjoint of the boundary operator, i.e. 
the coboundary operator, 8. The coboundary operator is the discrete representation 
of the topological derivatives grad, curl and div. Since dd = 0, it follows from a 
discrete Newton-Leibniz, Stokes and Gauss theorem that applying the coboundary 
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operator twice is always zero, 55 = (see J21 |9l). The coboundary operator also 
has matrix representations, G, C and D, that are the transpose of the connectivity 
matrices. We obtain the following topological sequence, 



^Ai'-^yAr, 



(8) 



where CG = and DC = 0. These matrices will explicitly appear in the final matrix 
system. An illustration of DC = is given in Figure [5] More details on the structure 
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Fig. 3 The action of twice the coboundary operator 5 on a vorticity d.o.f. has a zero net result on 
its surrounding volumes, because they all have both a positive and a negative contribution from its 
neighboring velocity faces. 



of geometry, orientation and degrees of freedom can be found in Gerritsma et al, 

0. 

Mimetic Operators. Let W = H{Q,curl), V = H{Q,div) and Q = L 2 (Q). The 
discretization of the flow variables involves a projection operator, %, from the com- 
plete vector spaces W, V and Q, to the discrete vector spaces W/,, V), and Q/,. Here the 
flow variables are expressed in terms of d.o.f. defined on fc-cells, and corresponding 
interpolation functions (also called basis-functions). The projection operator actu- 
ally consists of two steps, a reduction operator, $%, that integrates the flow variables 
on £-cells, and a reconstruction operator, J> , that interpolates the d.o.f. using the 
appropriate basis -functions. These mimetic operators were defined in [2, 10|. A 
composition of the two operators gives the projection operator 7£/, = J" o g$\ 

Reduction operator 2% is simply defined by integration. It possesses the following 
commutation relations, 



^grad = G^, ^curl = C^, ^div = D&. 



(9) 



The treatment of the reconstruction operator leaves some freedom, as long as it 
satisfies the following properties: be the right inverse of the reduction, &J? = Id, 
be the approximate left inverse of the reduction, J* 3?, = Id + tf(h p ), and it should 
possess the following commutation relations, 



gradJ^ = J^G, curlJ^ = J^C, div^^J^D. 



(10) 



1 For completeness, in a Hilbert setting the projection needs an additional smoothing argument. 
This step is ignored here to increase readibility. See 1 8 1 for more details. 
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When both the reduction and reconstruction operators commute with continu- 
ous and discrete differentiation, than also the projection operator Kt, possesses 
a commutation relation with differentiation. In case of the divergence operator, 
which is relevant to obtain a pointwise divergence-free solution, the commu- 
tation relation is given by, 



A\vK h = div J 3% = JVSt = J^div = K h A\v . 
The commutation relations in case of divergence are illustrated below, 



(11) 



V 



.y 



div 



* Q 



-> t. 
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v h 



-> r 



div 



> Qh- 



v 



V h 



div 



Q 



Xh 



div 



+ Qh- 



Since property ( [TT] ) also holds for the grad and curl, we obtain the following 
complex for discrete vector spaces, 



_ grad curl ., div _ 

®h — >W h — >V h — >Q h . 



(12) 



In practice we use J?V)S% from ( fTTj ) in computations. Relation (jTTJ implies 
among others that it satisfies the discrete LBB condition, 



j3/, := inf sup 

qh£Qhv h ev h 



(q h ,div\ h ) n 
Ikftllfillv/illv 



>j3>0, 



(13) 



where j3 is the inf-sup constant of the continuous problem, ([T|. Whereas the LBB 
condition is a measurement for numerical stability, the commutation relation indi- 
cates physical correctness of the numerical method. This last is a much stronger 
statement, which includes also the former. 

The conditions on the reconstruction operator have led to the construction of 
mimetic spectral element basis-functions, [5 , 10|. Since we use a tensor-based con- 
struction of point, line, surface and volume corresponding basis-functions, we only 
need nodal and edge interpolation functions. The nodal interpolation functions are 
the well-known Lagrange polynomials. The edge polynomials were derived from 
the Lagrange polynomials, based on the given conditions. For a set of Lagrange 
polynomials, k{x), i — Q,...,N, the edge polynomials, e,(x), i = 1 , . . . ,N, are given 

by, 



/.-=() 



6x 



(14) 



The Lagrange and edge polynomials possess the condition 8i<# = Id, i.e., 
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f x j 

li(xj) = 5ij, / ei(x)6x=Sij, (15) 

where Stj is the Kronecker delta. The interpolation function for a variable associated 
to a surface, for example, is given by, stj^ (x,y, z) = {k (x)ej (y)et (z) , e, (x)lj {y)ek (z) , 
ei(x)ej(y)k(z)}. 



Example 1 (Divergence operator in 2D). One of the most interesting proper- 
ties of the mimetic method presented in this paper, is that within our weak 
formulation, the divergence-free constraint is satisfied pointwise. Let m/, £ Vh 
be the velocity flux defined as 

Uh -\fiU4*wWj<y))- (6) 

Then the change of mass, m/, € Qh, is equal to the divergence of U/,, 

N N 

m h = div Uh = Y,Y, ( U 'J ~ u >- 1 J + V 'J ~ V 'J- 1 ) e i ( x ) e J b>) ■ 

N N 

where my = Ujj — Uj-ij + Vij — Vij-1 can be compactly written as m = Du. 
Note that if the mass production is zero, as in our model problem ( fTc] i, the 
incompressibility constraint can already be satisfied at the discrete level. In- 
terpolation using ei(x)ej(y) then results in a solution of velocity u/, that is 
pointwise divergence-free. 



5 A priori error estimates 

By standard interpolation theory it follows that we obtain the following /z-convergence 
rates for the interpolation errors of the flow variables, 

\\C0-7t h C0\\ H (cml) - &(h N ), ||u-%u|| ff(div) = 0(hP), \\p - 7t h p\\ L 2 = ff{h N ), 

(18) 
and that 1 1 div u — div %u|| L 2 = due to the commuting property. 

In cases with empty harmonic vector spaces, we have that the discrete vector 
spaces are conforming, i.e., W/, C W, V/ 7 C V and Qt, C Q. Moreover, due to the com- 
muting property, it follows that these spaces are compatible, i.e., curlW/j C V), and 
div Vi, = Qh- Finally they possess a Helmholtz-Hodge decomposition, a = grad0 + 
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curl*v and v = curia + grad*q'. In terms of vector spaces, this is, Wf, = Zw h ®Z^ 
and Vft = Zy h © Zy , where Z refers to the kernel or nullspace and Z x to its orthogo- 
nal complement. Having all these properties, a priori error estimates are derived in 
(8l that show optimal convergence rates for all admissible boundary conditions, in- 
cluding the no-slip boundary condition, which is non-trivial in mixed finite element 
methods. The a priori error estimates are given by 

\\(0-(O h \\ w <C inf IIo-ctaIIw, (19) 

o h eW h 

Mil — uJly < C inf llu — vJ|y+C inf \\co — Oh\\w> (20) 

Vh£V h a h GWi, 

\\p-Ph\\Q<C inf \\p-q b \\ Q + C inf ||u - \ h \\ v +C inf \\(D-G h \\ w , (21) 
q h eQ h v h eV h o h ew h 

where the constants C will differ in each case and are independent of h. It shows 
that the rate of convergence of the approximation errors are the same as those of the 
interpolation errors. 



6 Numerical Results 

For many years, the lid-driven cavity flow was considered one of the classical bench- 
mark cases for the assessment of numerical methods and the verification of incom- 
pressible (Navier)-Stokes codes. The 3D lid-driven cavity test case deals with a flow 
in a unit box with five solid boundaries and moving lid as the top boundary, mov- 
ing with constant velocity equal to minus one in x-direction. Especially the two line 
singularities make the lid-driven cavity problem a challenging test case. 

The left plot in Figure HI shows slices of the magnitude of the velocity field in a 
three dimensional lid-driven cavity Stokes problem, obtained on a 2 x 2 x 2 element 
mesh, where each element contains a Gauss-Lobatto mesh of Af = 8. The slices are 
taken at 10%, 50% and 90% of the y-axis. The right plot in Figure HJ shows slices of 
divergence of the velocity field. Figure HI confirms that the mixed mimetic spectral 
element method leads to an accurate result with a divergence-free solution. 

The second testcase shows the optimal convergence behavior for a 2D Stokes 
problem with no-slip boundary conditions. The testcase originates from a recent 
paper by Arnold et al [1], where sub-optimal convergence is shown and proven for 
no-slip boundary conditions when using Raviart-Thomas elements. Since Raviart- 
Thomas elements are the most popular //(div, £i) conforming elements, we compare 
our method to these results. 

Figure E\ shows the results of the Stokes problem on a unit square with velocity 
and pressure fields given by u = [ — 2x 2 (x—1 ) 2 y(2y — 1 ) (y — 1 ) , 2y 2 (y — 1 ) 2 x(2x — 

1 ) (x — 1 )] , p — (x— j) 5 + (y— 5 ) 5 ■ While for velocity both methods show optimal 
convergence, for pressure a difference of I is noticed in the rate of convergence and 
for vorticity a difference in rate of convergence of | is revealed. 
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Fig. 4 Left: slices of magnitude of the velocity field of a three dimensional lid-driven cavity Stokes 
problem obtained on a 2 x 2 x 2 element mesh with N = 8. Right: slices of the divergence of 
velocity. Is confirms a divergence-free velocity field. 



■ Raviart-Thomas, N=2 

■ Mimetic Spectral, N=2 
Mimetic Spectral, N=3 




Fig. 5 Comparison of the /i-convergence between Raviart-Thomas and Mimetic spectral element 
projections for the 2D Stokes problem with no-slip boundary conditions. 
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